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Building upon an analytical technique introduced by Chung and Peschel, Phys. Rev. B 64, art. 
064412 (2001), we calculated the many-body density matrix ps of a finite block of B sites within an 
infinite system of free spinless fermions in arbitrary dimensions. In terms of the block Green function 
matrix G (whose elements are Gy = ( c I c j)' wnere c \ an d c- are fermion creation and annihilation 
operators acting on sites i and j within the block respectively), the density matrix can be written as 
Pb = det(l — G) exp[£\ .(log G(t — G)~ x )ijc\c^\. Our results suggests that Hilbert space truncation 
schemes should retain the states created by a subset of the c|'s (in any combination), rather than 
selecting eigenvectors of pB independently based on the eigenvalue. 
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I. INTRODUCTION 

Exact solutions are hard to come by in many-body 
problems, and every so often we have to resort to nu- 
merical solutions. The traditional approaches, applied 
to finite systems, are exact diagonalization and quantum 
Monte Carlo. For quantum lattice models of fermions, 
the former is constrained by the size of the Hilbert space, 
which grows exponentially with the number of sites, while 
the latter is plagued by the 'minus-sign problem'. For 
quantum lattice models of bosons, the Hilbert space is 
infinite-dimensional even for finite systems. In either 
case, because of the enormous computational complexity 
involved, there is no hope of getting to the thermody- 
namic limit of infinite system size. In view of these diffi- 
culties, one then hopes for the next best thing: approxi- 
mate solutions that captures the essence of the physics. 

This is where renormalization group (RG) approaches 
comes in. In such approaches 1-11 to the approximate 
solution of otherwise intractable problems, the size of 
the Hilbert space is kept in check by aggressive trun- 
cation, with the hope that the small number of states 
kept will reproduce the more important features of the 
physics. Whatever the RG scheme, ultimately its success 
will lie in how the truncation is done. Since the quantum- 
mechanical state of a block embedded in a larger system 
must in general be described by a density matrix, it is 
therefore natural to use it to guide the truncations. 

With the goal of understanding the structure and spec- 
trum of the density matrix, and their implications on 
RG truncation schemes in mind, Peschel et al calcu- 
lated exactly the half-chain density matrix for several 
models. 12-14 For a chain of coupled harmonic oscillators 
and spinless Bogoliubov fermions, they found that the 
half-chain density matrices can be expressed exactly as 
the exponential of a pseudo-Hamiltonian, whose spec- 
trum is generated by a set of independent bosonic and 
fermionic operators respectively. In this paper, we pursue 



their analysis further for a system of free spinless fermions 
to obtain a closed-form formula relating the density ma- 
trix of a subsystem and the subsystem Green function 
matrix G (to be defined in Section II C). 

The organization of the paper will be as follows: we 
will start in Section II by reviewing the density matrix 
formulation of quantum mechanics, and how the density 
matrix ps of a subsystem can be obtained from the den- 
sity matrix po of the overall system. Following this, we 
will describe an alternative approach to calculating the 
density matrix elements as expectations of referencing 
operators. We shall show that the real-space structure, 
and the strong signs that point to a closed-form expres- 
sion for ps, is most readily discerned within this alter- 
native formulation. Then, in Section III, we derive this 
closed- form expression for ps in terms of the subsystem 
Green function matrix G by adapting the technique put 
forward by Chung and Peschel. 14 The existence of such a 
relation between pb and G tells us that ps is completely 
determined by its 0- and 1-particle sectors. We discuss 
the implications of this in Section IV, where we illustrate 
how the eigenvectors and eigenvalues of the (F > 1)- 
sectors of ps can be constructed out of the eigenvectors 
and eigenvalues of the 1-particle sector. We also show 
how symmetries of the Hamiltonian that are realized in 
ps affect the pattern of degeneracies in the eigenvalues 
of these sectors, an understanding of which is important 
in formulating a consistent truncation scheme. 



II. DENSITY MATRIX OF A FINITE BLOCK 

In this section, we first review the density matrix no- 
tions that will be used throughout this paper. Following 
this we develop the first of our two paths to calculate the 
density matrix ps for a particular block within a large 
system of non-interacting fermions. By analyzing the in- 
dex structure of the matrices involved, we arrived at a 
conjecture for a closed-form expression for the 1-particle 
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sector of the block density matrix ps in terms of the 
block Green function matrix G. 



A. Density Matrix Formulation of Quantum 
Mechanics 

In quantum mechanics one distinguishes between pure 
states, which occur, for example, at T — when the sys- 
tem is totally decoupled from the rest of the universe, 
and mixed states, which occur, for example, at T > 
when the system is in thermodynamic equilibrium with 
the rest of the universe. A pure state can be described by 
a wavefunction \^>) in the usual formulation of quantum 
mechanics, whereas a mixed state cannot. Both type of 
states are treated on equal footing in the density matrix 
formulation of quantum mechanics, in which the state of 
a system is described by a density matrix p (see, for ex- 
ample, Rcf. 15). In this formulation of quantum mechan- 
ics, the expectation of an operator A in a state described 
by po is given by 

(A) = Tt(p A). (2.1) 

If the state so described is pure, i.e. given in the usual 
formulation by the wavefunction |*), so that (A) = 
then it is clear that p = |*) (*|. 
In this paper, we shall be mainly interested in a finite 
subsystem of B sites, which we call the block, embedded 
within a larger system of N sites, with periodic boundary 
conditions in d dimensions. The latter can then be taken 
to the thermodynamic limit of infinite number of sites, 
i.e. N — > oo. The system minus the block is called the 
environment of the block. If the overall system is known 
to be in a pure state then in general the quantum- 
mechanical state of the block cannot be described by a 
pure state wavefunction. Instead, the mixed state of the 
block must be described by a block density matrix ps 
(see arguments in Ref. 16), so defined that 

(A) = Tt(p b A), (2.2) 

if the operator A acts entirely within the block. 

There are two useful formulas to relate the block den- 
sity matrix ps to the density matrix po of the entire sys- 
tem. The first formula, which we will used in Section III, 
follows from (2.1) and (2.2). Using the subscripts B and 
E respectively to make the trace over the degrees of free- 
dom associated with the block and its environment more 
explicit, we can rewrite (2.1) as 

(A) = Tt(p A) = Ti BiE {p A). (2.3) 

Since A does not act on the environment, we can trivially 
trace over environmental degrees of freedom to get 

(A)=Tr B {[Tr E (p )}A}. (2.4) 

Comparing this with (2.2), we find a consistent expecta- 
tion for A whether it is taken over the entire system or 



just over the block, if the block density matrix is defined 
as 

PB = Tr B (p ). (2.5) 

The second formula for ps allows us to write down its 
matrix elements explicitly when the overall system is in 
a pure state. To arrive at this formula, let us first note 
that any pure state of the overall system can be written 
as = J2b |k) l efc )> wnere l&) is a complete orthonormal 
(many-body) basis for the block, and |ej,) is the (unnor- 
malized) state of the environment associated with the 
state \b) on the block. Using this form for in (2.1), 
we find that 

(A) = J2 (eb\ (b\ A \b') \e v ) = Tt b { P bA) (2.6) 

6,6' 

if the block density matrix p B is defined such that 

(pB)b>b = (e b \ev) , (2.7) 

i.e. the matrix element of ps between \b) and \b') is none 
other than the overlap between their associated environ- 
mental states \eb) and \ey). 

B. Free Spinless Fermions 

Let us now apply (2.7) to calculate the block density 
matrix from the ground state of a ring of N — > oo free 
spinless fermions, the simplest realization of which is de- 
scribed by a translationally-invariant Hamiltonian with 
nearest-neighbor hopping 

ff = -*£( c h- +c 5 c i)' (2-8) 

<«> 

where c i and cj are the fermion annihilation and creation 
operators acting on site i, and {ij} runs once over each 
pair of neighbor sites. 

The Hamiltonian given in (2.8) is diagonal in momen- 
tum space, and can be written as 

# = E e 4<V (2-9) 

k 

Here 

J (2.10) 

c£ = JV-V^ c t e -ifcT, 

3 

are the momentum space annihilation and creation oper- 
ators, Ti is the position of site i, and Ck the single-particle 
energy associated with wavevector k. The ground state 
of such a system is just a Fermi sea 

i*f>= n 4io>, (2.H) 

k filled 
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where |0) is the vacuum, and the product is over the 
wavevectors inside the Fermi surface. 

As noted in (2.7), when the ground state wavefunction 
is written as \^f) — J2b 1^) \ e b)i the block density ma- 
trix elements are (pB)b'b = (eb\eb')- When dealing with 
a finite block and an infinite environment, it makes no 
sense to evaluate these environmental overlaps by first 
calculating |e&) and lev). Instead, we find that it possi- 
ble to evaluate such environmental overlaps with the help 
of operator products that are defined entirely within the 
block. To do so, let us first write the many-body states 
\b) on the block in the occupation number representation 



as | b) 



-,b„b 



where 



or 1 depending on 



whether the site j on the block is empty or occupied in 
the state We then define the referencing operators 



K b = Y\[n)c ] + {l-n))c 3 c]], 



(2.12) 



such that the effect of K b acting on a state \b') is Kb \b') = 
Sbb' |0) B , where \0) B is the reference state for which all 
sites on the block are empty. Letting Kb act on \^f) 
gives K b |*f> = Eb<< K b \b") \e b ») = \0) B \e b ). Hence, in 



the 



terms of the operators Kb and their conjugates K.\, 
density matrix elements are found to be 

( PB )b'b = (e b \e v ) = (9 F \KlK v \iS F ) = {K\K v ) . 

(2.13) 

From the way the operators Kb are defined, we know 
that p B is real and symmetric. Furthermore, (pB)b'b van- 
ishes if the states \b) and |6') do not contain the same 
number of fcrmions F. Consequently, the non-zero ma- 
trix elements of p B are found in a total of (B + 1) subma- 
trices along the diagonal, corresponding to the various 
-F-particle sectors, for F = 0, 1, . . . , B. We shall call 
such submatrices pb,f, and their eigenvalues the density 
matrix weights Wb,f,i, where 1 = 1,... , rank(pB ; j?). 



C. Conjecture Based on Index Structure 

In general, for a block of B sites, there are a total 
of 2 B Kb operators we need to write down explicitly to 
calculate the ~2 S density matrix elements. For large 
blocks, this is extremely tedious and has to be automated 
(see Appendix A), but for small blocks, it is not difficult 
to work out exact expressions for (pB)b'b in terms of the 
2n-point functions 



G 7 



(4c,) 



2(2-1) 



Gijklmn = ( c i c j c k C l C m c n) 



(-1) : 



C- C- C- 
C- C- C- 
Gjci Gj, m Gg ra 



(2.14a) 
(2.14b) 

(2.14c) 



and so forth, where i,j,k,l,--- = 1, . . . , B are sites on the 
block. As shown explicitly above, the 2n-point functions 
Gi 1 —i n j 1 —j n Wick factorizes into sums of products of 2- 
point functions Gjj for our non-interacting system, with 
an overall fermion factor of (— 1)"(™ -1 )/ 2 . 

At this point let us note that since the 2-point func- 
tions Gjj are labelled by two indices, it is convenient to 
organize them into a system Green function matrix ^ 
given by 



Gbi 

G~b + ii 



G N1 



Gib 

&2B 



G 
G 



1B+1 
2B+1 



G 



G 



BB 



G 



G 



BB+l 



B + 1B 



G NB 



^NB + 1 



Gin 

&2N 



G 



G- 



BN 



B+1N 



G NN 



of which 



G 



Gn G 12 ■ ■ ■ G 1B 
G21 G22 ■ ■ ■ G2B 



Gbi G B 2 • • • G 



BB 



(2.15) 



(2.16) 



is its restriction to the block. We call G the block Green 
function matrix. As a result of the translational invari- 
ance of H, & is also translationally invariant. In real 
space, this means that its matrix elements 5% = = 
(cjcj) are functions only of rj — rj. When ^ is restricted 
to the block to give G, however, this translational invari- 
ance is lost due to the fact that the presence of a block 
in the system allows an unambiguous definition of the 
origin. 

Anyway, from (2.12) and (2.13), we see on the one 
hand that (pB)bb' can be written as sums of 2n-point 
functions — which themselves factor into sums of prod- 
ucts of 2-point functions — and so we find that (pB)bb' 
are all functions of Gj, ; . On the other hand, the 1-particle 
sector of pb contains matrix elements (pB)bb' connecting 
the states \b) and \b'), which contain one particle each 
at sites, say, i and j respectively. Therefore, the matrix 
elements within p B \ may be indexed using i and j in- 
stead of b and b' . Diligently writing down the polynomial 
expressions 



{PB,i)ij = a irMh G kih 



fei, 



^ a ij^k 1 k2hhGk 1 l 1 Gk 2 l 2 + 



fei,fc 2 , 



^ a ijl 1 ---k B i 1 --i B Gk 1 i 1 x • • • x G 



k 1 ,...,k B , 
Ii,...,Ib 



k B l B ' 



(2.17) 



4 



we find that: (a) the coefficients ct[j) ki ... k ^...j are inde- 
pendent of i and j; and (b) indices other than i and j 
always appear in pairs, as if they are summed over. 

Exhaustively comparing the matrix elements of pB,\ 
and powers of G for 2 < B < 5, we find that 



P bs = G + G 2 -GTr(G) + 

G 3 - G 2 Tr(G) - | {Tr(G 2 ) - 
G 4 -G 3 Tr(G)-±{Tr(G 2 )- 
|iTr(G 3 ) - | Tr(G) Tr(G 2 ) 



[Tr(G)] 2 }G + 

[Tr(G)] 2 }G 2 - 

-i[Tr(G) 3 ]}G+-.. 

(2.18) 



What is most fascinating about this series is that for 
B = 2, (2.12) and (2.13) tell us that ps,\ can be at most 
0(G 2 ), since its matrix elements never contain terms 
with more than two creation and annihilation operators 
each. Yet (2.18) is perfectly valid for B = 2, because 
terms higher order in G vanish. For B = 3 and B = 4, 
we find similarly that terms higher order than 0(G 3 ) 
and 0(G 4 ) vanishes, respectively. If we conjecture that 

(2.18) gives the leading terms to an infinite series that 
holds true for all B > 5, then we can factorize it into 

Pb,i = (G + G 2 + G 3 + ■■■) x 

exp [-Tr(G+±G 2 + ±G 3 + •••)] • (2.19) 

Noting that the series inside the trace is just — log(l — G), 

(2.19) can be compactly written as 



p B ,\ = G(l-G)- 1 dct(l-G). 



(2.20) 



III. DERIVATION AND PROPERTIES OF qb 

In passing from (2.18) to (2.20), a leap of faith was re- 
quired, and it would appear forbiddingly difficult to actu- 
ally prove (2.20) for arbitrary block sizes B, by the alge- 
braic manipulations used in Section II C. Fortunately, an 
alternate technique introduced by Chung and Peschel 14 
can be adapted and extended for calculating the density 
matrix of a finite block, although it comes with its own 
set of technical difficulties. It turns out that if the whole 
system were in the Fermi sea ground state, the derivation 
would require the inversion of singular matrices. In the 
end, the singularities do cancel and give a well-defined 
answer, but a regularization is needed to avoid diver- 
gences in the intermediate steps. The most natural way 
to do so would be to generalize our problem to nonzero 
temperature, in which case the limit T — > then provides 
the needed regularization. 19 

In essence, the calculations is just that of evaluating a 
Gaussian integral with the usual shift in integration vari- 
ables. However, because we are dealing with fermions, 
whose creation and annihilation operators anticommute 
rather than commute, additional machinery is needed 
to accomplish the feat of Gaussian integration. After 
casting the system density matrix po as a Gaussian of 



the fermion operators, we introduce fermionic coherent 
states with the aid of anticommuting Grassmann vari- 
ables. The matrix elements of po between such coherent 
states, obtained via a translation machinery, are simi- 
larly of Gaussian form, but are now easier to handle. A 
Gaussian integration over the environmental degrees of 
freedom then yields elements of the block density matrix 
Pb, following which reverse translation gives ps proper. 



A. Exponential Form for po 

To get the calculations underway, we consider the 
grand-canonical T > density matrix p of the over- 
all system that the block is embedded in. As always, this 
is given by 



Po 



g- 1 exp[-[3(H - pF)} 



(3.1) 



where (i = l/fc^T, p is the chemical potential, and F = 
Sfe ^k^k ~ Si c \ c i i s tne fermion number operator. The 
prefactor i? -1 in (3.1) is just the reciprocal of the grand 
partition function, to ensure that Tr(po) = 1- 

The notations can be made more compact if we intro- 
duce the matrices V and its Fourier transform T, such 
that 

po = exp Fij4 c j) = ^ ex P (Efe r/cfe4 5 fc) . 

(3.2) 

where we have made use of the fact that H — p,F, and 
hence T, is diagonal in momentum space. The matrix 
elements of T can be read off from (2.8) as 



r 



(3p, if i = j; 

(it, if i and j are nearest neighbors; (3.3) 
0, otherwise, 



while those of T can be read off from (2.9) as 

ffcfc - -0E k , (3.4) 

where = — p is the single-particle energy measured 
relative to p. 

In order to prove our conjecture (2.20), it is clear that 
we need to somehow relate T to <£ '. To do this, let us 
note that since is translationally invariant, its Fourier 
transform 'S is diagonal in momentum space, with matrix 
elements given in the grand-canonical ensemble as 



%k = (c{c k ) 
observing which we find that 



1 



expf3E k + 1' 



'kk 



= exp(f fefc ) exp(f fcfe ) + 1 



(3.5) 



(3.6) 



But since both Sf and f are diagonal matrices, we have 
the relation 



e r = #(l-#)- 1 , 



(3.7) 
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where e r is the matrix exponential of T. 

Of course, and Sf corresponds merely to the matrix of 
the same Hilbert space operator evaluated in two diferent 
bases, and the same is true of T and V. As such, the 
matrix relation (3.7) between e r and holds true for e r 
and <S as well, i.e. we have 



(3.8) 



where £ = {£i , . . . , £b} are Grassmann coordinates asso- 
ciated with sites on the block, and r) = {771, . . . ,ry7v-s} 
are Grassmann coordinates associated with sites in the 
environment. 



Matrix Block Form 



B. Key Formulas Involving Grassmann Variables 

In the next stage of our derivations, we need to make 
use of Grassmann variables. These are anticommuting 
c-numbers familiar in the context of field theory (see for 
example, Ref. 17). If £j and £j are Grassmann variables, 
where i ^ j, then we have = — and £? = = 
£|. The purpose of introducing these is to define the 
fcrmionic coherent states 



IO = l£i6---6vHex P l-YZitel 



|0), (3.9) 



which are eigenstates of the fcrmion annihilation opera- 
tors, i.e. |£) = — & The value of coherent states in 
general is that one can replace the manipulation of non- 
commuting operators by the manipulation of c-number 
matrix elements. In the present case of fermions, anti- 
commuting operators may be made to commute by the 
insertion of Grassmann coefficients. 

There are three key formulas involving Grassmann al- 
gebra that we need for the derivations in this section. The 
first involves the matrix element of an exponentiated bi- 
linear operator exp . r^-ctc- j between fermionic co- 
herent states |£) and given by 



(£1 cxp (Eij r ijC Jc.) |0 = exp [E M (e r )^^ 

(3.10) 

where e r is the exponential of the matrix T. The second 
formula expresses the trace of an operator A as a Grass- 
mann integral over its coherent state matrix elements as 

Tr(A) = f Hd&dZie-X'^ (-£L4|£> . (3.11) 

The third formula that we would need is the Gaussian 
integral over Grassmann variables, 



•T.j, k £j A ik£k _ det A. 



(3.12) 



Seeing that po is written in (3.2) as the exponential of 
a quadratic form with coefficient matrix T, we make use 
of (3.10) to write down its matrix element between the 
fermionic coherent states |£t/) and \£' rj') as a Gaussian 
in Grassmann variables: 



(^IpoKVH^T'exp 



. (3.14) 



Our task now is to derive the matrix elements of ps in the 
same Gaussian form, after tracing out the environment. 

To find the matrix elements (£|/0b|£ ) of the density 
matrix ps on the block of B sites, we use (3.11) and 
perform a partial trace over the environment to give 

{Z\PB\g) = J drfdrie-*' 1 " (€ -r,\(>o\€ v) 
= £~ x j drj*dr) x 



exp 
exp 



v -+) (S !) ft 
(r -1*) e r a 



1 e"** 1 *' y d V * 



exp 



^77 x 



(3.15) 



Following this we must express these matrix elements in 
a form that would allow us to trace over the environment. 
To do so, let us first write (1 + e r ) in matrix block form 
as 



l + e r = 



A B 
B T C 



(3.16) 



The strategy then would be to evaluate the matrix el- 
ements of po in (3.2) using (3.10), follow the prescription 
in (2.5) where we trace over the environmental degrees of 
freedom using (3.11), and then use (3.10) in reverse to re- 
cover pb from its coherent state matrix elements. Before 
we do so, let us first tidy up the notations by relabelling 
the coherent states as 



\£v) = \£,i---£,b;vi---vn-b) 
= exp (- J2f=i & c 



£.7=1 ^ C 



10} , 



(3.13) 



where A is the B x B block submatrix, obtained by re- 
stricting the indices i and j of (1 + e r ) in coordinate 
space to range only over sites on the block, C is the 
(N — B) x (N — B) environment submatrix, obtained by 
restricting the indices i and j of (1 + e r ) to range only 
over sites in the environment, and B is the B x [N — B) 
decoherence submatrix of (1 + e r ), obtained by restrict- 
ing the the row index to range only over sites on the 
block and the column index to range only over sites in 
the environment. 
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D. Tracing Down p 

With (3.15) and (3.16), the block density matrix ele- 
ments can then be written as 

(Z\Pb\0 =^- 1 e«*( A - 1 )€' J dr t *dr t et* Br i-*>* BT t'- r i* c *>. 

(3.17) 

Here we have made use of the fact that since the Grass- 
mann variables occur quadratically in each term in the 
exponential, they commute with one another and we may 
thus factor the exponential as if it is an exponential of 
c-numbers. 

By performing a shift of the integration variables rj 
and 77*, and then evaluating the Grassmann Gaussian 
integral using (3.12), we find that 

(€IPb|€'> = £- 1 dctCet^ A - 1 - BC ~ lBT X', (3.18) 

which parallels Eq. (14) in Rcf. 14. From (3.18), we see 
that the expression for (£|ab|£') involves only the Grass- 
man coordinates £j and ^ associated with sites on the 
block. This is good. But it also involve the decoherence 
submatrix B as well as the environment submatrix G, 
with the latter appearing both in the exponential as well 
as in the normalization constant. These matrices have 
indices that range over sites outside the block, which we 
are supposed to have traced out and gotten over with. 

Indeed, this must have been successfully done, since 
A — 1 — BC~ 1 B T is a B x B matrix whose indices range 
only over sites on the block. In fact, using (B5) in Ap- 
pendix B, we can express this matrix entirely in terms of 
submatrices on the block, and write (3.18) as 



($\pB\$') = ^- 1 detCe^ D - 1 - 1 ^', 



(3.19) 



where D is the block submatrix of (1 + e r ) _1 , obtained 
by restricting its indices to range only over sites on the 
block. That leaves only the det G in the normalization 
that we have to deal with. 

To express =S _1 detG in terms of submatrices whose 
indices range only over sites on the block, we make use 
of the fact that 

Tr(p B ) = 1 = | dCdte-f* {-t\p B \Z) 

= J2~ 1 detC J dedSe-S' - 1 * ( 3 - 2 °) 
= .S^detCdet IT 1 , 



which means that 



^"MetG = dctL>. 



(3.21) 



With this we have succeeded in writing down a Gaus- 
sian form for the coherent state matrix elements of ps 
involving only degrees of freedom on the block. Using 



the translation machinery provided by (3.10), we then 
establish the exponentiated form 

PB = detD exp{]T [logOD- 1 - 1)] .. c\ Cj } (3.22) 



of Chung and Peschel. 

At this point, let us remark that the above formula 
for pb is of no practical use, if to find the matrix D, 
we actually have to evaluate the matrix (1 + e r ), whose 
indices run over the entire system, take its inverse (1 + 
e r ) _1 , and then from this identify the block submatrix 
D. This is essentially what was done in Rcf. 14, where 
the matrix parallel to D^ 1 — t was computed numerically, 
for the case of an environment equal in size to the block. 
For our problem, identifying A — BC~ 1 B T as D^ 1 with 
the aid of our analytic relation (B5) allows us to work 
with arbitrary, even infinite, environment sizes. 

Furthermore, armed with the relationship (3.8) ob- 
tained in Section III A, we can find that the normal- 
ization and matrix of coefficients appearing in (3.19) in 
terms of the block Green function matrix G. From 



1 + e r = 1 + Sf (1 - ST)" 1 = (1 - <S)- 



(3.23) 



we see that D is just (1 — Sf ) restricted to the block, 
i.e. D = 1-G, and consequently, D^ 1 = (l-G) -1 . With 
this, the normalization constant for ps can be written 
as det I? = det(l — G). For the matrix of coefficients 
(D^ 1 — 1) in the exponential, we see that 

D- 1 - 1 = (1 - G)" 1 - 1 = G(l - G)" 1 . (3.24) 

With this substitution, the matrix elements of ps now 
reads as 

(£\pb\£') = det(l - G) cxp [C G(l - G) -1 £'] . (3.25) 

so that, after using (3.10) in reverse translation, we can 
read off the operator form of ps as 



p B = det(l - G) cxp \ J2 [logG(l - G)- 1 ]^ c\ Cj 

3 

(3.26) 

In a suitable basis diagonalizing logG(l — G) _1 , this be- 
comes 



p B = det(l - G) cxp [- J2i Vififi 



(3.27) 



where the //'s are linear combinations of c^'s, and <pi is 
the associated pseudo-energy (see (3.33) for definition). 
With (3.27), we see that to find ps, we need only cal- 
culate the B x B block Green function matrix G from 
the ground state wavefunction with the aid of operators 
local to the block, and diagonalize it to determine fi and 
subsequently ipi. 

To connect this with the results that we obtained in 
Section II C, let us evaluate the matrix elements for the 
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0- and 1-particle sectors of ps- Taylor expanding the 
exponential in (3.27) gives us 

PB = det(l - G) J] [l + (e"*" - l)fln] , (3.28) 
i 

and so we see that the 0-particle sector is given by 

ps,o = B<0|p B |0) s =det(l-G), (3.29) 

while in the basis diagonalizing ps 1 the matrix elements 
in the 1-particle sector are given by 

s(0|/ ; ps,i//|0) B = det(l - G) B (0|/ ; l//|0) s + 

£(e- w -l)B<0|/i/ I t/f/ I t |0) J 
/' 

= det(l-G) e- Vl 

= det(l-G) [G{t-G)- 1 ] u . 

(3.30) 

This completes the proof of our conjecture at the end of 
Section II that as a matrix, ps,i is related to G by (2.20). 

E. The Pseudo-Energies ipi 

With the closed- form formula (3.27) for ps at hand, we 
arc now ready to understand its structure and spectra. 
To begin with, we find that the exponential form 



p B = det(l - G)exp 



= det(l - G)exp 



£(logG(l -Gr% c : Cj 



(3.31) 

where we define the matrix $ to be 

$ = -logG(l-G)" 1 = -logG + log(l-G), (3.32) 

implies that the weights and eigenvectors of the (F > 
l)-particle sectors of ps are determined completely by 
those in the 0- and 1-particle sectors. Defining the set of 
pseudo-energies 

tpi = -(logG(l-G)- 1 ) H , (3.33) 

for I = 1,...,B, which are the eigenvalues of <&, and 
ipa = — log dot (1 — G), we find that the weights of the 
1-particle block can be written as 

wi = cxp[-(^ + <pi)], (3.34) 
and pb can be written in the form 



ps = e Vo exp 



e- Vo exp[-H}. (3.35) 



Inspired by the resemblance of the form of ps to the 
density matrix of a quantum canonical ensemble, we call 
H the pseudo-Hamiltonian. 



F. Particle-Hole Symmetry at Half-Filling 

To complete our understanding of the structure and 
spectrum of pb, we want to know how symmetries of the 
original problem are built into ps- In particular, we will 
consider particle-hole symmetry on a bipartite lattice, 
on which we define a 'charge- conjugation' operator 
with "if 2 = I. 20 The action of ^ on the coordinate space 
fermion operators can be defined to be 



■to?- 



(3.36) 



where ( — 1)* is defined to be +1 (—1) whenever the 
site i belongs to the even (odd) sublattice. In a d- 
dimcnsional hypercubic lattice, where the site index is 

i = {h,i2,---, "id}, the factor (— l) 1 is rightfully given by 

(— 1)* = (_i)«i+*2H — H<i. 

There are two conditions, one on the dispersion re- 
lation efc, and the other on the chemical potential p,, 
implied by particle-hole symmetry. To derive the first 
condition on the dispersion relation, we note from (3.36) 
that in momentum space — when the lattice is a Bravais 
lattice — that 



-k+Q> 



(3.37) 



where the wavevector Q is defined by e l( 3' ri = (— 1)\ 21 
We can then check, with (2.9) and (3.37), that 



k+Q^-k-Q 



^ e -k'-QCk'+2Q^k' 



Now, from the definition of Q, it is clear that 



C k'+2Q 



,i(k'+2Q)- rj 



= N- 1 ^^2e ik '^[(-l) j ] 2 c j 



(3.38) 



(3.39) 



and thus (dropping the prime on the dummy wavevector 
k' that is summed over) 

W = e -*-Q~ c kA = e - k -Q + £ - € -k-Qc\ck- 

k k k 

(3.40) 

For time-reversal invariant systems, e_fc = efe. Also, 
for our choice of Hamiltonian, J2k e -k-Q = J2k' e *=' = 
TriJ = 0. Thus (3.40) simplifies to 



tf/fif = £-e fc+0 4c fc . 



(3.41) 



Since it is clear from (2.8) and (3.36) that tfHtf = H, 
(3.41) tells us that the dispersion relation associated with 
the particle-hole symmetric Hamiltonian H must satisfy 
the condition 



£fe+Q 



(3.42) 



Next, to understand how the second condition on the 
chemical potential comes about, let us note the trivial 
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fact that, since is a reduced density matrix of po, 
for there to be any sense in talking about the manifes- 
tation of particle- hole symmetry in ps, po must first be 
particle- hole symmetric, i.e. ^ po& = po- When p is the 
density matrix of the ground state at T = 0, then it is 
particle- hole symmetric whenever the ground state \^f) 
is. For I^f) to be particle-hole symmetric, it must have 
the transformation property 

^N=WN, (3.43) 

where r/<g — ±1 is a phase factor associated with ^ . We 
know that this is satisfied only by the half-filled ground 
state. At finite temperature, when po is taken from the 
grand canonical ensemble and has the form given in (3.1), 
what, if any, extra conditions must be satisfied in order 
for po to be particle-hole symmetric? 

Indeed, there appears to be cause for concern: un- 
like H, which is invariant under 'charge-conjugation', the 
fermion number operator F transforms under as 

«W = Vc\ Ci V = Y,Ci4 = N- F, (3.44) 

i i 

and so for p to be particle-hole symmetric, i.e. 

tfpo^ = exp 0[H - p{N - F)] = p 0l (3.45) 

we must have p, = 0. For a dispersion relation satisfying 
(3.42), p = corresponds to precisely the situation of 
half-filling. At least for the grand canonical ensemble, 
there appears to be no other conditions necessary for p 
to be particle-hole symmetric. 

With these conditions in mind, we may now proceed to 
investigate how particle-hole symmetry shows up in the 
pseudo-energy spectrum (and hence the spectrum of the 
block density matrix ps)- But first, we must understand 
how the symmetry is manifested in the Green function 
matrix £f , and its restriction to the block, G. Knowing 
from our arguments above that p = 0, we see that the 
matrix elements of <S in momentum space simplifies to 

&kk = ^——r- (3.46) 

exp/3e fc + 1 

Furthermore, using (3.42), we can relate k+Q to &kk 

by 

^ +Q ' fe+Q = exp/3 £fc+Q + l = cxp(-/3 efe ) + 1 = 1 -^- 

(3.47) 

This gives rise to the condition 

S% - % - (-lf-^% 3 (3.48) 

that must be satisfied by the matrix elements of Sf in 
coordinate space. 

This same result can be derived more transparently for 
the special case of T = 0: using the fact that c tf 2 = 1, 



rj 2 * = 1, as well as (3.36) and (3.43), we find that 
S% = (* F | c t Cj .|* F > 

= (-iy+i {y F \ Cl c]\y F ) (3.49) 

= (-iy+i Sij - (-iy+j (* F \c] Ci \* F ) 

= Sij — (—1) ^ Jc $iji 

where we have made use of the fact that Sf is symmetric, 
i.e. — • 

Since (3.49) is a condition satisfied by the matrix el- 
ements of ^ individually, it holds just as well to those 
restricted to the block, i.e. Gij. In particular, we note 
that (3.49) can actually be written as a matrix equation, 
which when restricted to the block reads as 

G = 1 - JGJ, (3.50) 

where J = diag(e I( ^' ri ) = diag(l, — 1, 1, — 1, . . . ) in coor- 
dinate space, and J 2 = 1. 

To appreciate the implications of (3.50), let us consider 
an eigenvector |A/) of G correspond to the eigenvalue A;, 
such that 

G\\ l ) = \ l \\ l ). (3.51) 

By (3.32), this is also the eigenvector of pB, with corre- 
sponding pseudo-energy 

W = -logAi+log(l-Ai). (3.52) 
Using (3.50), we find that 

GJ \X t ) = (1 - JGJ)J \X{) 
= J\\ l )-JGJ 2 \\ l ) 
= J\\i)-JG\\i) (3.53) 
= J|A;)-A;J|A*) 
= (1-Ai)J|Ai), 

i.e. the state |A;/) = J|A/) generated by particle-hole 
symmetry from |A;) is also an eigenvector of G, with 
eigenvalue A// = (1 — A ; ). The pseudo-energy tpy as- 
sociated with | A;/) is then 

<pi> = - log A;/ + log(l - Xu) = -(pi. (3.54) 

It is interesting to note how the weights Wb,i,i, being ex- 
ponentials of the pseudo-energies <pi , hide this particular 
aspect of particle-hole symmetry. 

IV. THE (F > 1)-PARTICLE SECTORS 

Up to this point, our discussions have been for ar- 
bitrary dimensions. To demonstrate how the (F > 1)- 
particle sectors can be constructed from the 0- and 1- 
particle sectors, we specialize to the 1-dimensional case, 
wherein the Fermi sea is 

\k\=n{-K/a) 

\*f)= J! 410), (4-1) 

|A;|=0 
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where a is the lattice constant and n is the filling fraction. 
The 2-point functions can be computed explicitly as 



Gu - 



simrn\i — j\ 
n\i-j\ 



(4.2) 



We now illustrate how to construct the weights and eigen- 
vectors of the (F > l)-particle sectors of ps starting from 
(po, ipi and using the example of B = 3 at half-filling. 

At half-filling, n = \, the 2-point functions Gjj take 
on particularly simple values 



Gn — G 77 — G 33 — 75, G 



2 ' 



'12 — G 2 3 — \, G 13 — 0, 
(4-3) 

with which we find, using the machinery developed in 
Section II B, the 0-particle and 1-particle sectors of p 3 to 
be 



P3,o - (000|p 3 |000) = § - 

(100|p 3 |100) (100|p 3 |010) <100|p 3 |001)- 

(010|p 3 |100) (OlOlpglOlO) (010|p 3 |001) 

(001|p 3 |100) (001|p 3 |010) (001|p 3 |001). 



P3,l 



(4.4) 



r 1 

8 




1 

In 


1 -1 

7? 


1 


1 




1 


2w 


8 




2lT 


1 




1 

2ir 


1 

8 J 



Diagonalizing these, we find 

™3,o,i = (^-i) (75 + £)' 



"■<■ 17s- - 



W3,l,2 



\V8 



VV8 



W 3 

and thus 

Pa = 



(4.5a) 



(4.5b) 



<Pi = - log 



— + - 
_i 1 



<p 2 = ~ log 1 = 0, 



¥>3 = - log 



1 

78 



78 + 



-2.94777. . 



+2.94777. . 



(4.6a) 



(4.6b) 



Since ipi = — <p 3 , we call these a particle-hole conjugate 
pair of pseudo-energies, and say that P3 is the particle- 
hole conjugate of p\. The eigenvectors of the 1-particle 
sector are 



t 



fl 



fl = + 



i r t , j_ r t , i t 

2 C 1 ' ^2 C 2 * 2 C 3> 
_ r t _ J_ r t 



4 



v/2 C 2 



(4.7a) 
(4.7b) 
(4.7c) 



corresponding to </?i, <^2 and tp$ respectively. 

We can easily check that the fl's obey the same an- 
ticommutator relation as the cj's, i.e. they obey Pauli's 
Exclusion Principle, and so the eigenvectors of the 2- 
particle sector of p 3 are created by 



flfl 



2 C 2 C 3 



V2 



t„t lJJ 



ci c 



1°3 



2 c 1 c 2 , 



J1J3 - ^2 2 3 v^ C l C 2' 



C 2 C 3 + 7/9 1 3 2 L 1 C 2 



(4.8a) 
(4.8b) 
(4.8c) 



with associated pseudo-energies p\ + P2 = Pi, Pi + p>3 = 
= (fi2 and tp2 + P3 = P3 respectively. Here we see 
that because of the particle-hole symmetry in the ground 
state wavefunction, the pseudo-energies of the 2-particle 
sector are identical to those in the 1-particle sector, which 
implies that the density matrix weights of the 2-particlc 
sector are identical to those in the 1-particle sector. 
For the 3-particle sector, we find that the eigenvec- 



tftft _ J J J 

-1 L 2 L 3' 



associ- 



tor is created by the operator flflf: 
ated with pseudo-energy ipi + ip 2 + p>3 = 0, and hence 
^3,3,1 = e ~ Vo = w 3,o.i- This method of generating all 
(F > l)-particle sectors, starting from the 0- and 1- 
particle sectors, for larger block sizes at various filling 
fractions n was verified numerically. 

Another manifestation of particle-hole symmetry is a 
queer degeneracy between weights in the F- and (F + 2)- 
particle sectors. This we understand as follows: if ipi ± + 
••• + ipi F is a weight in the -F-particle sector, then in gen- 
eral we can find weights (Lp h -\ \-<Pl F ) + <Pl F +i+<PlF+2 = 

{ipi t +• • - + (pi F ) in the (_F+2)-particle sector, where pi F+1 
and pi F+2 are particle-hole conjugates of each other. 

In fact, from the construction outlined above, we also 
know the pattern of degeneracy. For example, for B = 6 
at half-filling, the pseudo-energies are of the form — ip a , 
—ifb, —<Pc, fc <Pb and ip a , corresponding to the eigen- 
statcs created by /I , /l b , fl c , /], fl and /] respec- 
tively, where tp a > ipb > <p c - We then see in the 3-particle 
sector that flf^_ b f- a |0) and fcf- c f-a |0) have the same 
pseudo-energy of ip b - ip b - ip a = (p a = ip c - (p c + ip ai and 
are thus degenerate, whereas fcf-bf-a |0) I s nondegen- 
erate with pseudo-energy ip c — <f b — 'Pa- 



V. CONCLUSIONS 

To summarize, in this paper we showed that elements 
of the block density matrix, (pB)b'b, can be calculated 
as the expectation (K b Kb') of a product of referencing 
operators K b and Kb>, which are themselves formed out 
of fcrmion operators Cj and cj local to the block. By 
inspecting the matrix elements (pB,i)ij and Gjj of the 1- 
particle sector of ps and the block Green function matrix 
G respectively for block sizes up to B = 5, we are led to 
a conjecture of the closed-form relation (2.20) between 
Pbs and G. 

Adapting the technique that Chung and Peschel used 
to calculate the half-chain density matrix of a chain of 
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spinlcss Bogoliubov fcrmions, we find that we can not 
only prove this conjecture, but also derive a closed-form 
relation (3.26) between the entire block density matrix 
Pb and G, which can also be written in (3.27) as the 
exponential of a pseudo-Hamiltonian H . The spectrum 
of H is generated by the independent fermion operators 
f}, which also generate the eigenvectors of G, and hence 
can be determined by diagonalizing G. It is amusing 
to numerically compute the pseudo-Hamiltonian in real 
space. For n ^ \, the generic form of H (a bilinear in 
{cj} and {q}) admits hopping to all other sites on the 
block, as well as an on-site potential. The symmetry at 
half-filling ensures that the on-site potential is zero and 
hopping only connects to the other sublattice. 

We identify three important implications of (3.27) in 
formulating truncation schemes based on ps, for the pur- 
pose of performing an RG analysis. Firstly, we note 
that the spectrum of ps is completely determined by the 
block Green function matrix G. It suffice therefore to 
calculate the eigenvectors and eigenvalues in its 0- and 
1-particle sectors. The eigenvectors and eigenvalues of 
the (F > 1) -particle sectors of pb, if needed, can be 
systematically generated from the fermion operators fj 
and their pseudo-energies ipi , as illustrated in Section IV. 
This fact is evident in Rcf. 14, but its significance was not 
emphasized. If one is studying the density matrix of a 
noninteracting toy model (as in this paper), we have an 
enormous saving in terms of computational time: instead 
of diagonalizing the entire pb, which is of rank 0(e B ) 1 
we can diagonalize just the 1-particle sector pb,i, which 
is of rank 0(B). Possible objects of such a study could 
be: (i) the distribution of eigenvalues; 18 (ii) errors in the 
dispersion relation due to truncation; 18 and (iii) compar- 
ing the product basis of two blocks of length B with the 
basis of one block of length 2B, to weigh the effects of 
the correlations respectively neglected or included. 

Secondly, it is highly desirable in RG calculations to 
ensure that the truncation scheme preserves the symme- 
tries of the target state. Using the specific example of 
particle-hole symmetry, we saw in Section IV that pb.f 
and pb,b-f have the same set of weights, and the eigen- 
vectors of pb,f are related, up to a phase, to the eigenvec- 
tors of pb.b-f acted upon by So . Naively, we might ex- 
pect that to preserve particle-hole symmetry, all we have 
to do is to keep SS \wB,F.k) in pb,b-f if \wB,F,k) in Pb,f is 
kept. However, there is more to particle-hole symmetry. 
Under the action of So ', the half-filled ground state \^f) 
goes (up to a phase) back to itself. Within the block, this 
global symmetry transformation brings the mixed state 
of the block back to the same mixed state. Because the 
mixed state of the block does not have a definite particle 
number, particle-hole symmetry is not merely a relation 
between pb^f and pb,b-f- Rather, particle-hole sym- 
metry imposes strict conditions on the spectra of pb.f 
and pb,f', for F, F' = 0, . . . , B. In fact, in Section IIIF, 
we elaborated on the condition that particle-hole sym- 
metry imposes on the 1-particle sector. This condition 
is most intuitive when written in terms of the eigenval- 



ues A; of G or the pseudo-energies ipi, but not immedi- 
ately apparent if we just stare at the 1-particle density 
matrix weights wbs,i- It is therefore dangerous to base 
symmetry-preserving truncation schemes on ps and its 
eigenvalues alone. 

This brings us to the last of the implications that we 
wish to highlight. While a toy noninteracting model is 
studied in this paper, our ultimate goal is to address in- 
teracting systems, particularly Fermi liquids. Since these 
(in their low-energy limit) have the same eigcnstate struc- 
ture as a noninteracting Fermi sea (after a unitary trans- 
formation), their density matrices also should have the 
same structure as a noninteracting system. The explicit 
form of the many-body density matrix, as exhibited in 
Section III of this paper, hints at the proper design of 
truncation schemes. Rather than independently truncat- 
ing in each F'-particle sector, we should define the trun- 
cated states using a set of 'creation operators' which sat- 
isfy the usual anticommutation relations, and quite likely 
these are closely related to the approximate quasiparti- 
cle creation operators, which should be constructed as 
a product of the renormalization scheme. We will have 
more discussions on the implications of such a trunca- 
tion scheme based on picking out a set of appropriate 
'creation operators' for the numerical study of interact- 
ing systems, the role of dimensionality, and comparisons 
with the conventional DMRG, in a second paper. 18 

Based on our observations on the pattern of degen- 
eracies within and between the F-particle sectors of ps 
in Section IV, we realize that if the truncation is car- 
ried out naively, there is a very real danger of ending up 
with an inconsistent scheme of truncation. This prob- 
lem occurs quite generally, at various filling fractions and 
block sizes, but can be most clearly illustrated using our 
example of B = 6 at half-filling. For example, let us 
say that as the result of a naive truncation, the states 
f-a |0> , /1 6 |0), /l c |0) and fl |0) in the 1-particle sector 
are kept. Examining the 2-particle sector, we find the 
states fL a fa |0) and /l c /j |0), which are degenerate in 
their pseudo-energies. We can build up the latter, but 
not the former, using the 1-particle operators kept, and 
so we should keep the latter but not the former. If we 
truncate the 2-particle sector naively, then based on the 
density matrix weights alone we would be probably end 
up keeping or throwing out both of these 2-particle states. 

In fact, the situation for naive truncation is worse, 
since the state f_ a fl |0) has lower pseudo-energy than 

f-cfc |0) an d b e kept instead. We see therefore that 
naive truncation is likely to led to inconsistencies: some 
many-particle states built up from the 1-particle states 
kept get thrown out, while other many-particle states 
that cannot be built up from the set of 1-particle states 
kept end up being retained. Hence, we find that as far as 
particle-conserving models are concerned, for any trun- 
cation scheme to be consistent, the truncation must be 
carried out on the 1-particle sector of pb alone. 

Finally, let us remark that everything done in this 
paper can be trivially extended to the case of spinfull 
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fermions, so long as they are noninteracting. Every ob- 
ject in our calculations, in particular the Fermi sea wave- 
function and the block density matrix, will merely in the 
spinfull case be replaced by the direct product of two 
such objects with spin-up and spin-down flavors. 



The limiting consideration for the whole calculation 
is the diagonalization time, which is determined by the 
condition number of pb, rather than array storage. In 
general, the condition number, which is the ratio of the 
largest weight to the smallest west, grows exponentially 
with system size. 18 
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APPENDIX A: AUTOMATING THE 
NUMERICAL COMPUTATION OF p B 

As we saw in Sec. II, each of the 2 B basis configura- 
tions of the block corresponds to an operator Kb, so that 
(ps)b'b = {K^Kv). Therefore, to obtain all matrix ele- 
ments of it sufficed to automate the calculation of ex- 
pectations of an arbitrary string of creation/annihilation 
operators (each operator acting on one site). 

First, this formal string of operators must be simpli- 
fied. Through a systematic set of anticommutations, it is 
brought to a canonical form, such that (a) it is normal- 
ordered, with one substring of all creation operators fol- 
lowed by one substring of all annihilation operators; and 
(b) within each substring the operators are ordered by 
the site. Of course, each site can appear at most once in 
each substring (otherwise it reduces to zero.) A compli- 
cation of this step is that the result is generally a sum of 
many terms in the canonical form, since every rearrange- 
ment of the form c x c\ — > 1 — c{c 1 produces two terms 
from one. 

Next, we note that within this sum, only terms contain- 
ing a balanced number, say n, each of creation and anni- 
hilation operators will contribute to the expectation. By 
the Wick theorem, such 2n-point functions Gi 1 —i n j 1 —j n 
reduces to the determinant of an n x n matrix, as shown 
in (2.14). 

For models in which fermion number F is conserved, we 
can further separate ps into the various i^-particle sec- 
tors pb,f before diagonalization. This is particularly im- 
portantly at half-filling, for there exists generic degenera- 
cies between states in different sectors (see Section IV), 
and there is a danger that a naive diagonalization of the 
whole matrix ps will produce eigenstates with mixed par- 
ticle number. 



APPENDIX B: BLOCK INVERSION FORMULA 

Consider a square N x N symmetric matrix ^ written 
in matrix block form as 



A B 

B T C 



(Bl) 



where A is a square i\q x N\ symmetric matrix, B is a 
Ni x N 2 non-square matrix and C is a square N 2 x N 2 
symmetric matrix. Here Ni + N 2 = N. 

If we write the inverse matrix also in the matrix 

block form 



D E 
E T F 



(B2) 



where D is a square iVi x Ni symmetric matrix, E is a 
Ni x N 2 non-square matrix and F is a square N 2 x N 2 
symmetric matrix, how are D, E and F related to the 
matrix blocks A, B and C in ./#? 

Using the fact that Ji = 1, and thus 



(B3) 



(where the subscripts, which will henceforth be dropped 
for notational clarity, following the 1's and D's indicate 
the shape and size of the matrices) we find the following 
relations between the matrix blocks of Jt and Ji~ x : 



' A 


B 




D 


E 




IjVixATi 




B T 


C 




E T 


F 






1n 2 xn 2 _ 



AD + BE T = 1, 
AE + BF = (D, 
B T D + CE T = 0, 
B T E + CF = 1. 



(B4a) 
(B4b) 

(B4c) 
(B4d) 



Solving for D, E and F in terms of A, B and C, we find 
that 

D=[A-BC- 1 B T ]~ 1 , (B5a) 
E=-A- 1 B(C-B T A- 1 By 1 , (B5b) 
F = [C - B T A- 1 By 1 . (B5c) 
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n 



\k\<k F 



clckU\ k '\>k F £ k'cl- This is pos- 



sible only if Fkk = +oo for \k\ < fc^ and Tkk = — oo 
for | A; | > ItF- For the purpose of algebraic manipula- 
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as e = e 



A 1 + (e A - e A ) «f . It is then easy to show that 

A\-l /i i „-A\-ll</? 



(1 + e r )- 1 = (1 + e^)- 1 ! + [(1 + e A )- x - (1 + e 
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